Overview

This analysis examines the effects of wounding on respiration rates in two coral species: - Porites sp. (por) - Acropora pulchra (acr)

Experimental Design: - 36 total corals (18 per species) - 3 treatments per species (n=6 per treatment): - Control (0): No wounding - Wound Type 1 (1): First wound method - Wound Type 2 (2): Second wound method - Wound date: 2023-05-27 - Measurements over 4 timepoints spanning 23 days

Physiological Measurements: 1. Respiration rates (O2 consumption) 2. Photosynthesis rates (O2 production) 3. Growth (buoyant weight) 4. Photosynthetic efficiency (Fv/Fm via PAM fluorometry) 5. Surface area (wax dipping method)

Load Libraries and Functions

suppressPackageStartupMessages({
  library(tidyverse)
  library(lmerTest)
  library(emmeans)
  library(rstatix)
  library(janitor)
  library(ggthemes)
  library(knitr)
  library(kableExtra)
  library(pbkrtest)
  library(sjPlot)
  library(cowplot)
  library(patchwork)
  library(here)
})

# Source custom functions (loads remaining packages without conflicts)
source("scripts/analysis_functions.R")

# Set theme for all plots
theme_set(theme_few(base_size = 12))

# Define consistent color scheme and labels for all plots
treatment_colors <- c("0" = "#2E86AB", "1" = "#A23B72", "2" = "#F18F01")
treatment_labels <- c("0" = "Control", "1" = "Small Wound", "2" = "Large Wound")

# Helper functions that may need Rmisc
if (!requireNamespace("Rmisc", quietly = TRUE)) {
  install.packages("Rmisc")
}

Load and Prepare Data

Sample Information

# Load sample metadata
sample_info <- read_csv("data/metadata/sample_info.csv", show_col_types = FALSE)

# Convert to factors
sample_info <- sample_info %>%
  mutate(
    treatment = as.factor(treatment),
    genus = as.factor(genus),
    wound_date = as.Date(as.character(wound_date), format = "%Y%m%d")
  )

# Summary table
sample_info %>%
  group_by(genus, treatment) %>%
  summarize(n = n(), .groups = "drop") %>%
  pivot_wider(names_from = treatment, values_from = n, names_prefix = "Treatment_") %>%
  kable(caption = "Sample sizes by genus and treatment") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Sample sizes by genus and treatment
genus	Treatment_0	Treatment_1	Treatment_2
acr	6	6	6
por	6	6	6

Growth Data

# Load processed weight data
initial_weight <- read_csv("data/processed/growth/initial_weight.csv", show_col_types = FALSE)
postwound_weight <- read_csv("data/processed/growth/postwound_weight.csv", show_col_types = FALSE)
day7_weight <- read_csv("data/processed/growth/day7postwound_weight.csv", show_col_types = FALSE)
final_weight <- read_csv("data/processed/growth/final_weight.csv", show_col_types = FALSE)

# Combine all weights
growth_data <- initial_weight %>%
  select(coral_id, species, initial) %>%
  left_join(postwound_weight %>% select(coral_id, species, postwound),
            by = c("coral_id", "species")) %>%
  left_join(day7_weight %>% select(coral_id, species, day7),
            by = c("coral_id", "species")) %>%
  left_join(final_weight %>% select(coral_id, species, final),
            by = c("coral_id", "species"))

# Calculate growth rates
growth_data <- growth_data %>%
  mutate(
    growth_g = final - postwound,
    growth_rate_g_day = growth_g / 23,  # 23 days post-wounding
    growth_normalized = growth_rate_g_day / initial
  )

# Add treatment info
growth_data <- growth_data %>%
  left_join(sample_info %>% select(coral_id, genus, treatment),
            by = c("coral_id", "species" = "genus"))

# Summary
growth_data %>%
  group_by(species, treatment) %>%
  summarize(
    mean_growth = mean(growth_rate_g_day, na.rm = TRUE),
    se_growth = sd(growth_rate_g_day, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  ) %>%
  kable(caption = "Growth rates by species and treatment (g/day)", digits = 4) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Growth rates by species and treatment (g/day)
species	treatment	mean_growth	se_growth
acr	0	0.0743	0.0066
acr	1	0.0511	0.0087
acr	2	0.0516	0.0108
por	0	0.1910	0.0208
por	1	0.0872	0.0665
por	2	0.0535	0.0192

PAM Data (Fv/Fm)

# Load averaged Fv/Fm data
pam_data <- read_csv("data/processed/pam/average_fvfm.csv", show_col_types = FALSE) %>%
  select(-any_of("X")) %>%  # Drop X column if it exists
  mutate(
    treatment = as.factor(treatment),
    genus = as.factor(genus),
    date = as.factor(date)
  )

# Convert date to timepoint label
pam_data <- pam_data %>%
  mutate(
    timepoint = case_when(
      date == "20230603" ~ "Day 7",
      date == "20230619" ~ "Day 23",
      TRUE ~ as.character(date)
    ),
    timepoint = factor(timepoint, levels = c("Day 7", "Day 23"))
  )

# Summary
pam_data %>%
  group_by(genus, treatment, timepoint) %>%
  summarize(
    mean_fvfm = mean(average_fv_fm, na.rm = TRUE),
    se_fvfm = sd(average_fv_fm, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  ) %>%
  kable(caption = "Fv/Fm by genus, treatment, and timepoint", digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Fv/Fm by genus, treatment, and timepoint
genus	treatment	timepoint	mean_fvfm	se_fvfm
acr	0	Day 7	697.222	10.119
acr	0	Day 23	705.500	9.017
acr	1	Day 7	702.000	8.572
acr	1	Day 23	704.278	6.034
acr	2	Day 7	699.333	9.146
acr	2	Day 23	710.167	5.626
por	0	Day 7	664.833	5.221
por	0	Day 23	648.611	11.821
por	1	Day 7	676.556	10.345
por	1	Day 23	652.500	7.172
por	2	Day 7	658.611	5.269
por	2	Day 23	656.778	5.349

Surface Area Data

# Load surface area data
sa_data <- read_csv("data/processed/surface_area/final_surface_areas.csv", show_col_types = FALSE) %>%
  select(-any_of(c("X", "...1"))) %>%
  rename(coral_id = coral_number, genus_sa = taxa) %>%
  mutate(genus_sa = tolower(genus_sa))

# Summary by genus
sa_data %>%
  group_by(genus_sa) %>%
  summarize(
    mean_sa = mean(CSA_cm2, na.rm = TRUE),
    sd_sa = sd(CSA_cm2, na.rm = TRUE),
    min_sa = min(CSA_cm2, na.rm = TRUE),
    max_sa = max(CSA_cm2, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  kable(caption = "Surface area summary (cm²)", digits = 2) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Surface area summary (cm²)
genus_sa	mean_sa	sd_sa	min_sa	max_sa
acropora	17.88	5.10	11.55	30.27
porites	24.49	6.98	15.28	43.76

Respirometry Data

# Load respiration rate data for Porites
# Note: Acropora data needs to be processed first

# Function to load rates from a date
load_rates <- function(date, rate_type, genus = "Porites") {
  path <- paste0("data/processed/respirometry/", genus, "/", rate_type, "/", date, "/")

  if (!dir.exists(path)) {
    return(NULL)
  }

  csv_files <- list.files(path, pattern = "\\.csv$", full.names = TRUE)

  if (length(csv_files) == 0) {
    return(NULL)
  }

  data <- read_csv(csv_files[1], show_col_types = FALSE) %>%
    mutate(date = date, rate_type = rate_type) %>%
    mutate(coral_id = as.numeric(coral_id))  # Standardize coral_id type

  return(data)
}

# Load all Porites respiration data
resp_dates <- c("20230526", "20230528", "20230603", "20230619")

resp_data_list <- list()
for (date in resp_dates) {
  resp <- load_rates(date, "Respiration", "Porites")
  photo <- load_rates(date, "Photosynthesis", "Porites")

  if (!is.null(resp)) resp_data_list[[paste0(date, "_resp")]] <- resp
  if (!is.null(photo)) resp_data_list[[paste0(date, "_photo")]] <- photo
}

# Combine all data
if (length(resp_data_list) > 0) {
  resp_all <- bind_rows(resp_data_list) %>%
    clean_names() %>%
    select(-starts_with("x")) %>%
    filter(!is.na(coral_id))

  # Convert date to timepoint
  resp_all <- resp_all %>%
    mutate(
      timepoint = case_when(
        date == "20230526" ~ "Pre-wound",
        date == "20230528" ~ "Day 1",
        date == "20230603" ~ "Day 7",
        date == "20230619" ~ "Day 23",
        TRUE ~ date
      ),
      timepoint_num = case_when(
        date == "20230526" ~ 0,
        date == "20230528" ~ 1,
        date == "20230603" ~ 7,
        date == "20230619" ~ 23,
        TRUE ~ as.numeric(NA)
      ),
      timepoint = factor(timepoint, levels = c("Pre-wound", "Day 1", "Day 7", "Day 23"))
    )

  # Add genus and treatment info
  resp_all <- resp_all %>%
    mutate(genus = "por") %>%
    left_join(sample_info %>% select(coral_id, treatment),
              by = "coral_id")

  cat("Loaded", nrow(resp_all), "respirometry measurements for Porites\n")
  cat("Dates:", paste(unique(resp_all$date), collapse = ", "), "\n")
  cat("Rate types:", paste(unique(resp_all$rate_type), collapse = ", "), "\n")
} else {
  cat("No respirometry data found\n")
  resp_all <- NULL
}

## Loaded 304 respirometry measurements for Porites
## Dates: 20230526, 20230528, 20230603, 20230619 
## Rate types: Respiration, Photosynthesis

Exploratory Data Analysis

Growth Patterns

# Growth by treatment and species
p1 <- ggplot(growth_data, aes(x = treatment, y = growth_rate_g_day, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5, size = 2) +
  facet_wrap(~species, labeller = labeller(species = c("acr" = "Acropora", "por" = "Porites"))) +
  scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
  labs(
    title = "Growth Rate by Treatment and Species",
    x = "Treatment",
    y = "Growth Rate (g/day)"
  ) +
  theme_few(base_size = 12)

# Normalized growth
p2 <- ggplot(growth_data, aes(x = treatment, y = growth_normalized, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5, size = 2) +
  facet_wrap(~species, labeller = labeller(species = c("acr" = "Acropora", "por" = "Porites"))) +
  scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
  labs(
    title = "Size-Normalized Growth Rate",
    x = "Treatment",
    y = "Normalized Growth Rate (g/day / initial mass)"
  ) +
  theme_few(base_size = 12)

print(p1)

print(p2)

# Save plots
ggsave("reports/Figures/growth_by_treatment.png", p1, width = 10, height = 6, dpi = 300)
ggsave("reports/Figures/growth_normalized.png", p2, width = 10, height = 6, dpi = 300)

PAM Fluorometry (Fv/Fm)

# Fv/Fm over time by treatment
p3 <- ggplot(pam_data, aes(x = timepoint, y = average_fv_fm, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5, size = 2) +
  facet_wrap(~genus, labeller = labeller(genus = c("acr" = "Acropora", "por" = "Porites"))) +
  scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
  labs(
    title = "Photosynthetic Efficiency (Fv/Fm) Over Time",
    x = "Timepoint",
    y = "Fv/Fm"
  ) +
  theme_few(base_size = 12)

# By treatment across species
p4 <- ggplot(pam_data, aes(x = genus, y = average_fv_fm, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  facet_wrap(~timepoint) +
  scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
  labs(
    title = "Fv/Fm by Species and Treatment",
    x = "Genus",
    y = "Fv/Fm"
  ) +
  scale_x_discrete(labels = c("acr" = "Acropora", "por" = "Porites")) +
  theme_few(base_size = 12)

print(p3)

print(p4)

ggsave("reports/Figures/fvfm_timeseries.png", p3, width = 10, height = 6, dpi = 300)
ggsave("reports/Figures/fvfm_by_species.png", p4, width = 10, height = 6, dpi = 300)

Respirometry Data

if (!is.null(resp_all)) {
  # Filter out blank corals (coral_id == 0 or 1)
  resp_filtered <- resp_all %>%
    filter(!coral_id %in% c(0, 1))

  # Respiration over time
  resp_only <- resp_filtered %>% filter(rate_type == "Respiration")

  p5 <- ggplot(resp_only, aes(x = timepoint, y = umol_l_sec, fill = treatment)) +
    geom_boxplot(alpha = 0.7) +
    geom_jitter(width = 0.2, alpha = 0.3, size = 1) +
    scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
    labs(
      title = "Respiration Rate Over Time (Porites)",
      subtitle = "Negative values indicate O2 consumption",
      x = "Timepoint",
      y = "Respiration Rate (µmol O2 / L / sec)"
    ) +
    theme_few(base_size = 12)

  # Photosynthesis over time
  photo_only <- resp_filtered %>% filter(rate_type == "Photosynthesis")

  p6 <- ggplot(photo_only, aes(x = timepoint, y = umol_l_sec, fill = treatment)) +
    geom_boxplot(alpha = 0.7) +
    geom_jitter(width = 0.2, alpha = 0.3, size = 1) +
    scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
    labs(
      title = "Photosynthesis Rate Over Time (Porites)",
      subtitle = "Positive values indicate O2 production",
      x = "Timepoint",
      y = "Photosynthesis Rate (µmol O2 / L / sec)"
    ) +
    theme_few(base_size = 12)

  print(p5)
  print(p6)

  ggsave("reports/Figures/respiration_timeseries.png", p5, width = 10, height = 6, dpi = 300)
  ggsave("reports/Figures/photosynthesis_timeseries.png", p6, width = 10, height = 6, dpi = 300)

  # Combined P and R plot
  p7 <- ggplot(resp_filtered, aes(x = timepoint_num, y = umol_l_sec, color = treatment)) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.5, position = position_dodge(width = 1)) +
    stat_summary(fun = "mean", geom = "point", size = 3, position = position_dodge(width = 1)) +
    stat_summary(fun = "mean", geom = "line", position = position_dodge(width = 1), aes(group = treatment)) +
    facet_wrap(~rate_type, scales = "free_y") +
    scale_color_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
    labs(
      title = "O2 Exchange Rates Over Time (Porites)",
      x = "Days Post-Wounding",
      y = "Rate (µmol O2 / L / sec)"
    ) +
    theme_few(base_size = 12)

  print(p7)
  ggsave("reports/Figures/O2_rates_combined.png", p7, width = 12, height = 6, dpi = 300)
}

Statistical Analysis

Growth Analysis

# Separate analysis for each species
cat("\n=== PORITES GROWTH ANALYSIS ===\n")

## 
## === PORITES GROWTH ANALYSIS ===

por_growth <- growth_data %>% filter(species == "por")

# Linear model
mod_growth_por <- lm(growth_normalized ~ treatment, data = por_growth)
print(anova(mod_growth_por))

## Analysis of Variance Table
## 
## Response: growth_normalized
##           Df     Sum Sq    Mean Sq F value Pr(>F)  
## treatment  2 0.00010679 5.3393e-05  3.8337 0.0452 *
## Residuals 15 0.00020891 1.3927e-05                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(summary(mod_growth_por))

## 
## Call:
## lm(formula = growth_normalized ~ treatment, data = por_growth)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0027700 -0.0019869 -0.0007704  0.0002741  0.0120689 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.008133   0.001524   5.338 8.28e-05 ***
## treatment1  -0.004687   0.002155  -2.175   0.0460 *  
## treatment2  -0.005540   0.002155  -2.571   0.0213 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.003732 on 15 degrees of freedom
## Multiple R-squared:  0.3383, Adjusted R-squared:   0.25 
## F-statistic: 3.834 on 2 and 15 DF,  p-value: 0.0452

# Pairwise comparisons
emm_growth_por <- emmeans(mod_growth_por, ~treatment)
pairs_growth_por <- pairs(emm_growth_por, adjust = "tukey")
print(pairs_growth_por)

##  contrast                estimate      SE df t.ratio p.value
##  treatment0 - treatment1 0.004687 0.00215 15   2.175  0.1082
##  treatment0 - treatment2 0.005540 0.00215 15   2.571  0.0525
##  treatment1 - treatment2 0.000853 0.00215 15   0.396  0.9177
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

cat("\n=== ACROPORA GROWTH ANALYSIS ===\n")

## 
## === ACROPORA GROWTH ANALYSIS ===

acr_growth <- growth_data %>% filter(species == "acr")

mod_growth_acr <- lm(growth_normalized ~ treatment, data = acr_growth)
print(anova(mod_growth_acr))

## Analysis of Variance Table
## 
## Response: growth_normalized
##           Df     Sum Sq    Mean Sq F value  Pr(>F)  
## treatment  2 6.0109e-05 3.0054e-05  4.0222 0.03994 *
## Residuals 15 1.1208e-04 7.4722e-06                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(summary(mod_growth_acr))

## 
## Call:
## lm(formula = growth_normalized ~ treatment, data = acr_growth)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0039210 -0.0025156 -0.0004061  0.0023428  0.0036627 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.010643   0.001116   9.537 9.29e-08 ***
## treatment1  -0.003605   0.001578  -2.285   0.0373 *  
## treatment2  -0.004100   0.001578  -2.598   0.0202 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002734 on 15 degrees of freedom
## Multiple R-squared:  0.3491, Adjusted R-squared:  0.2623 
## F-statistic: 4.022 on 2 and 15 DF,  p-value: 0.03994

emm_growth_acr <- emmeans(mod_growth_acr, ~treatment)
pairs_growth_acr <- pairs(emm_growth_acr, adjust = "tukey")
print(pairs_growth_acr)

##  contrast                estimate      SE df t.ratio p.value
##  treatment0 - treatment1 0.003605 0.00158 15   2.285  0.0891
##  treatment0 - treatment2 0.004100 0.00158 15   2.598  0.0500
##  treatment1 - treatment2 0.000495 0.00158 15   0.313  0.9475
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

PAM Analysis

# Linear mixed model with timepoint as fixed effect and coral as random effect
cat("\n=== PORITES FV/FM ANALYSIS ===\n")

## 
## === PORITES FV/FM ANALYSIS ===

por_pam <- pam_data %>% filter(genus == "por")

mod_pam_por <- lmer(average_fv_fm ~ treatment * timepoint + (1|coral_id), data = por_pam)
print(anova(mod_pam_por, type = 3))

## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## treatment            206.01  103.00     2    15  0.4192 0.66503  
## timepoint           1773.35 1773.35     1    15  7.2175 0.01691 *
## treatment:timepoint  762.23  381.11     2    15  1.5511 0.24417  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Pairwise comparisons
emm_pam_por <- emmeans(mod_pam_por, ~treatment|timepoint)
print(emm_pam_por)

## timepoint = Day 7:
##  treatment emmean   SE   df lower.CL upper.CL
##  0            665 7.98 26.6      648      681
##  1            677 7.98 26.6      660      693
##  2            659 7.98 26.6      642      675
## 
## timepoint = Day 23:
##  treatment emmean   SE   df lower.CL upper.CL
##  0            649 7.98 26.6      632      665
##  1            652 7.98 26.6      636      669
##  2            657 7.98 26.6      640      673
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

cat("\n=== ACROPORA FV/FM ANALYSIS ===\n")

## 
## === ACROPORA FV/FM ANALYSIS ===

acr_pam <- pam_data %>% filter(genus == "acr")

mod_pam_acr <- lmer(average_fv_fm ~ treatment * timepoint + (1|coral_id), data = acr_pam)
print(anova(mod_pam_acr, type = 3))

## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment            34.45   17.23     2    15  0.0632 0.9390
## timepoint           457.48  457.48     1    15  1.6792 0.2146
## treatment:timepoint 115.73   57.86     2    15  0.2124 0.8110

emm_pam_acr <- emmeans(mod_pam_acr, ~treatment|timepoint)
print(emm_pam_acr)

## timepoint = Day 7:
##  treatment emmean   SE df lower.CL upper.CL
##  0            697 8.26 27      680      714
##  1            702 8.26 27      685      719
##  2            699 8.26 27      682      716
## 
## timepoint = Day 23:
##  treatment emmean   SE df lower.CL upper.CL
##  0            706 8.26 27      689      722
##  1            704 8.26 27      687      721
##  2            710 8.26 27      693      727
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Respirometry Analysis

if (!is.null(resp_all)) {
  # Filter Porites, remove blanks, remove pre-wound timepoint
  resp_analysis <- resp_filtered %>%
    filter(timepoint != "Pre-wound")

  # Separate analysis for Respiration
  cat("\n=== PORITES RESPIRATION ANALYSIS ===\n")
  resp_resp <- resp_analysis %>% filter(rate_type == "Respiration")

  # Linear mixed model
  mod_resp <- lmer(umol_l_sec ~ treatment * timepoint_num + (1|coral_id),
                   data = resp_resp,
                   contrasts = list(treatment = "contr.sum"))
  print(anova(mod_resp, type = 3))

  # Pairwise comparisons by timepoint
  cat("\nPairwise comparisons at each timepoint:\n")
  for (tp in unique(resp_resp$timepoint)) {
    cat("\n", tp, ":\n")
    resp_tp <- resp_resp %>% filter(timepoint == tp)
    if (nrow(resp_tp) > 0) {
      mod_tp <- lm(umol_l_sec ~ treatment, data = resp_tp)
      emm_tp <- emmeans(mod_tp, ~treatment)
      print(pairs(emm_tp, adjust = "tukey"))
    }
  }

  # Analysis for Photosynthesis
  cat("\n=== PORITES PHOTOSYNTHESIS ANALYSIS ===\n")
  resp_photo <- resp_analysis %>% filter(rate_type == "Photosynthesis")

  mod_photo <- lmer(umol_l_sec ~ treatment * timepoint_num + (1|coral_id),
                    data = resp_photo,
                    contrasts = list(treatment = "contr.sum"))
  print(anova(mod_photo, type = 3))

  # Pairwise comparisons
  cat("\nPairwise comparisons at each timepoint:\n")
  for (tp in unique(resp_photo$timepoint)) {
    cat("\n", tp, ":\n")
    photo_tp <- resp_photo %>% filter(timepoint == tp)
    if (nrow(photo_tp) > 0) {
      mod_tp <- lm(umol_l_sec ~ treatment, data = photo_tp)
      emm_tp <- emmeans(mod_tp, ~treatment)
      print(pairs(emm_tp, adjust = "tukey"))
    }
  }
}

## 
## === PORITES RESPIRATION ANALYSIS ===
## Type III Analysis of Variance Table with Satterthwaite's method
##                             Sum Sq    Mean Sq NumDF  DenDF F value Pr(>F)
## treatment               1.5530e-05 7.7651e-06     2 97.386  0.7807 0.4609
## timepoint_num           1.7065e-05 1.7065e-05     1 85.905  1.7157 0.1937
## treatment:timepoint_num 9.9933e-06 4.9967e-06     2 85.905  0.5024 0.6069
## 
## Pairwise comparisons at each timepoint:
## 
##  Day 1 :
##  contrast                 estimate      SE df t.ratio p.value
##  treatment0 - treatment1 -0.000253 0.00037 33  -0.684  0.7742
##  treatment0 - treatment2  0.000245 0.00037 33   0.661  0.7877
##  treatment1 - treatment2  0.000498 0.00037 33   1.345  0.3811
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## 
##  Day 7 :
##  contrast                estimate      SE df t.ratio p.value
##  treatment0 - treatment1  0.00277 0.00223 33   1.245  0.4357
##  treatment0 - treatment2 -0.00121 0.00223 33  -0.542  0.8515
##  treatment1 - treatment2 -0.00398 0.00223 33  -1.787  0.1896
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## 
##  Day 23 :
##  contrast                 estimate       SE df t.ratio p.value
##  treatment0 - treatment1 -0.000617 0.000322 33  -1.918  0.1495
##  treatment0 - treatment2  0.000285 0.000322 33   0.886  0.6527
##  treatment1 - treatment2  0.000903 0.000322 33   2.804  0.0222
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## 
## === PORITES PHOTOSYNTHESIS ANALYSIS ===
## Type III Analysis of Variance Table with Satterthwaite's method
##                             Sum Sq    Mean Sq NumDF  DenDF F value Pr(>F)
## treatment               4.0666e-05 2.0333e-05     2 96.201  1.5586 0.2157
## timepoint_num           2.7186e-05 2.7186e-05     1 83.987  2.0839 0.1526
## treatment:timepoint_num 8.7910e-06 4.3957e-06     2 83.987  0.3369 0.7149
## 
## Pairwise comparisons at each timepoint:
## 
##  Day 1 :
##  contrast                 estimate       SE df t.ratio p.value
##  treatment0 - treatment1 -4.15e-05 0.000245 33  -0.169  0.9843
##  treatment0 - treatment2  2.83e-04 0.000245 33   1.156  0.4875
##  treatment1 - treatment2  3.25e-04 0.000245 33   1.325  0.3917
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## 
##  Day 7 :
##  contrast                 estimate     SE df t.ratio p.value
##  treatment0 - treatment1 -0.004650 0.0024 33  -1.934  0.1450
##  treatment0 - treatment2  0.000857 0.0024 33   0.357  0.9325
##  treatment1 - treatment2  0.005506 0.0024 33   2.291  0.0712
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## 
##  Day 23 :
##  contrast                 estimate       SE df t.ratio p.value
##  treatment0 - treatment1 -0.000330 0.000367 33  -0.901  0.6436
##  treatment0 - treatment2 -0.000168 0.000367 33  -0.458  0.8914
##  treatment1 - treatment2  0.000163 0.000367 33   0.443  0.8977
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Summary Tables

Sample Summary

# Overall sample summary
sample_info %>%
  group_by(genus, treatment) %>%
  summarize(n = n(), .groups = "drop") %>%
  pivot_wider(names_from = genus, values_from = n) %>%
  kable(caption = "Sample sizes by treatment and genus",
        col.names = c("Treatment", "Acropora", "Porites")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Sample sizes by treatment and genus
Treatment	Acropora	Porites
0	6	6
1	6	6
2	6	6

Mean Rates Summary

if (!is.null(resp_all)) {
  # Calculate mean rates by treatment, timepoint, and rate type
  summary_stats <- resp_filtered %>%
    group_by(treatment, timepoint, rate_type) %>%
    summarize(
      n = n(),
      mean_rate = mean(umol_l_sec, na.rm = TRUE),
      se_rate = sd(umol_l_sec, na.rm = TRUE) / sqrt(n()),
      .groups = "drop"
    )

  summary_stats %>%
    pivot_wider(
      names_from = rate_type,
      values_from = c(mean_rate, se_rate),
      names_sep = "_"
    ) %>%
    kable(caption = "Mean O2 exchange rates (µmol/L/sec) ± SE", digits = 4) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
    scroll_box(width = "100%", height = "400px")
}

Mean O2 exchange rates (µmol/L/sec) ± SE
treatment	timepoint	n	mean_rate_Photosynthesis	mean_rate_Respiration	se_rate_Photosynthesis	se_rate_Respiration
0	Pre-wound	12	-0.0011	-0.0031	0.0003	0.0005
0	Day 1	12	-0.0012	-0.0033	0.0002	0.0003
0	Day 7	12	0.0007	-0.0037	0.0005	0.0013
0	Day 23	12	0.0005	-0.0028	0.0002	0.0003
1	Pre-wound	12	-0.0001	-0.0029	0.0003	0.0002
1	Day 1	12	-0.0012	-0.0031	0.0002	0.0002
1	Day 7	12	0.0053	-0.0065	0.0029	0.0024
1	Day 23	12	0.0008	-0.0022	0.0002	0.0003
2	Pre-wound	12	0.0000	-0.0029	0.0003	0.0002
2	Day 1	12	-0.0015	-0.0036	0.0002	0.0003
2	Day 7	12	-0.0002	-0.0025	0.0003	0.0002
2	Day 23	12	0.0006	-0.0031	0.0003	0.0002

Conclusions

Key Findings

Growth

[Summary of growth results to be filled in after running analysis]

Photosynthetic Efficiency (Fv/Fm)

[Summary of PAM results to be filled in after running analysis]

Respiration

[Summary of respiration results to be filled in after running analysis]

Photosynthesis

[Summary of photosynthesis results to be filled in after running analysis]

Next Steps

Process Acropora respirometry data
Calculate normalized rates (per gram dry mass or per cm² surface area)
Calculate P:R ratios
Examine correlations between physiological measures
Compare temperature effects if variable temperatures present

Analysis completed: 2025-10-28 R version: R version 4.5.1 (2025-06-13)

Wound Effects on Respiration in Porites and Acropora pulchra

Regeneration Wound Respiration Project

2025-10-28

Overview

Load Libraries and Functions

Load and Prepare Data

Sample Information

Growth Data

PAM Data (Fv/Fm)

Surface Area Data

Respirometry Data

Exploratory Data Analysis

Growth Patterns

PAM Fluorometry (Fv/Fm)

Respirometry Data

Statistical Analysis

Growth Analysis

PAM Analysis

Respirometry Analysis

Summary Tables

Sample Summary

Mean Rates Summary

Conclusions

Key Findings

Growth

Photosynthetic Efficiency (Fv/Fm)

Respiration

Photosynthesis

Next Steps